home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
basic
/
chaosexe.zip
/
XINITFLO.TRU
< prev
next >
Wrap
Text File
|
1980-01-01
|
4KB
|
112 lines
!PROGRAM TITLE "XINITFLO"
CLEAR
PRINT" ***PENDULUM - INITIAL FLOW OF BLOCK OF COORDS***"
PRINT"THIS PROGRAM DISPLAYS THE EVOLUTION OF A BLOCK OF PHASE POINTS PENDULUM"
PRINT"AT VARIOUS TIMES DURING THE DRIVE CYCLE. FOR EXAMPLE, IF PHI=0.5 THEN"
PRINT"THE BLOCK WILL BE PRINTED EVERY 1/2 DRIVE PERIOD UNTIL THE MAXIMUM TIME."
PRINT"AS A FIRST TRIAL LET THE MINIMUM TIME BE ZERO AND THE MAXIMUM BE 5,"
PRINT" WHICH IS ABOUT 1/2 A DRIVE PERIOD."
LIBRARY "SGLIB.TRC"
!
DECLARE DEF accel
DIM A(1),B(1)
INPUT prompt"Input driving force strength: ":g
INPUT prompt"Input damping (If no damping then input 9999999):":q
!INPUT prompt"Input initial angle, angular velocity: ": xint,vint
INPUT prompt"input lower and upper xint:":lowerxint, upperxint
INPUT prompt"input lower and upper vint:":lowervint, uppervint
INPUT Prompt"Input min. and max. time:":tmin,tmax
INPUT prompt"Input phase angle/(2*pi): ":phi
CALL PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
CALL SETXSCALE(XMIN,XMAX)
CALL SETYSCALE(YMIN,YMAX)
CALL SETTEXT("PENDULUM INITIAL FLOW","ANGLE","ANGULAR VELOCITY")
CALL RESERVELEGEND
!
DATA O,O
CALL DATAGRAPH(A,B,1,O,"WHITE")
CALL gotocanvas
LET phi=phi*2*pi
FOR xint=lowerxint to upperxint step 0.03
FOR vint=lowervint to uppervint step 0.03
LET x=xint
LET v=vint
LET t=0
CALL graphpoint(xint,vint,1)
!
!CALCULATION AND GRAPHING BLOCK
FOR i=1 to 1000000
CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q) ! Take a step of size tstep
LET tshalf=tstep/2
CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q) !Take two half steps
CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
LET d1=abs(xn-xnew)
LET d2=abs(vn-vnew)
LET delta=max(d1,d2)
IF delta<eps then
IF t>tmin then
LET tnew=t+tstep
LET w1=mod(phi-w*t,2*pi) !Check for Poincare section
LET w2=mod(w*tnew-phi,2*pi)
IF w1<w*tstep then
IF w2<w*tstep then
LET ts=w1/w
CALL rk4(x,v,ts,xp,vp,t,w,g,q) !CALCULATES POINT AT SECTION
! IF abs(xp)>pi then LET xp=xp-2*pi*abs(xp)/xp
CALL GRAPHPOINT(XP,VP,1)
END IF
END IF
END IF
LET x=xnew
LET v=vnew
LET t=t+tstep !Expand step size
LET tstep=tstep*.95*(eps/delta)^.25
! IF abs(x)>pi then !bring theta back into range
! LET x=x-2*pi*abs(x)/x
! END IF
ELSE !else reduce step size
LET tstep=tstep*.95*(eps/delta)^.2
END IF
IF t>tmax then LET i=1000001
NEXT i
NEXT vint
NEXT xint
LET G$=STR$(G)
LET Q$=STR$(Q)
CALL ADDLEGEND("G="&STR$(G)&" Q="&STR$(Q)&" TIME="&STR$(PHI*3/2),0,1,"WHITE")
CALL DRAWLEGEND
get key variable
clear
END
!
!
SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
DECLARE DEF accel
LET xk1=tstep*v
LET vk1=tstep*accel(x,v,t,w,g,q)
LET xk2=tstep*(v+vk1/2)
LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
LET xk3=tstep*(v+vk2/2)
LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
LET xk4=tstep*(v+vk3)
LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
END SUB
DEF accel(x,v,t,w,g,q)
LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
END DEF
!
SUB PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
LET W=0.66666666
LET EPS=1.0E-6
LET TSTEP=0.5
LET XMIN=-6
LET XMAX=6
LET YMIN=-4
LET YMAX=4
END SUB